function ivp3a

%  solve IVP using MATLAB routines
%  y'' + sin(y)  with   y(0) = ep  y'(0)=0   '

clear *
clf
hold on


ep=1/3;
y0=ep;
tmax=200;

[t,y] = ode45(@rhs,[0 tmax],[y0 0]);

% get(gcf)
set(gcf,'Position', [1281 735 577 197]);

tt=linspace(0,tmax,800);
for it=1:800
	asy0(it)=cos(tt(it));
	asy1(it)=asy0(it)+ep*(tt(it)*sin(tt(it))/16-cos(3*tt(it))/192+cos(tt(it))/192);
	asy2(it)=cos(tt(it)*(1-ep^2/16));
end;

hold on
box on
yy=y(:,1)/ep;
plot(t,yy)
plot(tt,asy0,'--r','LineWidth',1)
%plot(tt,asy1,'-.k','LineWidth',1)
%plot(tt,asy2,'-.g','LineWidth',1)
%axis([0 2.5 0 0.601])
%title('\epsilon = 0.01','FontSize',14,'FontWeight','bold')
%legend(' Exact',' Asymptotic (1 term)','Location','South');
%legend(' Exact',' 1 term',' 2 terms','Location','EastOutside');
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off

figure
set(gcf,'Position', [1279 461 577 197]);
hold on
box on
yy=y(:,1)/ep;
plot(t,yy)
plot(tt,asy0,'--r','LineWidth',1)
%plot(tt,asy2,'-.k','LineWidth',1)
axis([180 200 -1 1])
%title('\epsilon = 0.01','FontSize',14,'FontWeight','bold')
%legend(' Exact',' 1 term',' 2 terms','Location','EastOutside');
legend(' Numerical',' Asymptotic (1 term)','Location','EastOutside');
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off


function dy=rhs(t,y)
dy=zeros(2,1);
ep=0.5;
dy(1)=y(2);
dy(2)=-sin(y(1));